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ABSTRACT 

> : 

Analytical studies and numerical simulations of time dependent axially 
■<^|- ! symmetric flows onto black holes have shown that it is possible to produce 

stationary shock waves with a stable position both for ideal inviscid and for 
00 ! moderately viscous accretion disks. 

We perform several two dimensional numerical simulations of accretion 
flows in the equatorial plane to study shock stability against non-axisymmetric 
azimuthal perturbations. We find a peculiar new result. A very small 
perturbation seems to produce an instability as it crosses the shock, but after 
some small oscillations, the shock wave suddenly transforms into an asymmetric 
j> \ closed pattern, and it stabilizes with a finite radial extent, despite the inflow 

^ j and outflow boundary conditions are perfectly symmetric. 

The main characteristics of the final flow are: 1) The deformed shock rotates 
steadily without any damping. It is a permanent feature and the thermal energy 
content and the emitted energy vary periodically with time. 2) This behavior 
is also stable against further perturbations. 3) The average shock is still very 
strong and well defined, and its average radial distance is somewhat larger than 
that of the original axially symmetric circular shock. 4) Shocks obtained with 
larger angular momentum exhibit more frequencies and beating phenomena. 5) 
The oscillations occur in a wide range of parameters, so this new effect may 
have relevant observational consequences, like (quasi) periodic oscillations, for 
the accretion of matter onto black holes. Typical time scales for the periods are 
0.01 and 1000 seconds for black holes with 10 and 10 6 solar mass, respectively. 
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1. Introduction 

Shock waves in rotating accretion flows onto compact objects present a very important 
mechanism to transform the potential gravitational energy into radiation. The possibility 
that they exist also around black holes (BH) has been suggested years ago by Hawley, 
Smarr, & Wilson ( 1984a, 1984b). However their shocks were unstable in radial direction. 
Recently it has been shown that there is a wide range of parameters (namely specific 
angular momentum and initial temperature) for which steady shock configurations are 
possible (Chakrabarti 1989, 1990; for viscous flows see Chakrabarti & Molteni 1993, 
Lanzafame et al. 1998). Essentially the shock position depends on thermal energy and 
angular momentum contents of the flow: the balance between the centrifugal and total 
pressure forces (ram plus internal) is the main factor determining the shock formation. The 
stability was examined with axially symmetric numerical calculations (Molteni, Lanzafame, 
& Chakrabarti 1994; Chakrabarti & Molteni 1995). Recent analytical studies (Lu & Yuang 
1997; Nakayama 1993, 1994) have also shown the stability of isothermal and adiabatic 
shocks against axially symmetric perturbations. 

The stability of shock waves in accretion flows has great astrophysical relevance since, 
in this way, the emission due to the shock is not a transient episode, but can be a permanent 
mechanism responsible of the radiated energy. However, the studies made up to now are all 
in axial symmetry. 

In this paper we examine the azimuthal stability of the initially axisymmetric shocks 
by perturbing them with an asymmetric perturbation. We explore the parameter space 
defined by the angular momentum and temperature of the accreted gas. The computations 
are done in the equatorial plane assuming slab symmetry in the z direction (thin disk 
approximation). We use two different computer codes to integrate numerically the 
hydrodynamic equations: the Versatile Advection Code (VAC, Toth 1996, 1997, also see 
http: / /www. phys.uu.nl/ "toth/p and another one developed by Kuznetsov (Kuznetsov et. 



al. 1998). 

Here we are not discussing the origin of the supersonic accretion flow and the influence 
of viscosity. The formation of supersonic flows in the accretion disk has been predicted by 
various authors in different contexts valid both for galactic and extragalactic supermassive 
Black Holes (for viscous transonic flows cfr. Chakrabarti 1996, for advection dominated 
flows see Gammie and Popham 1998, Igumenshev, Abramowicz, & Novikov 1998, Narayan, 
Mahadevan, & Quataert 1998). These works clearly show that since the viscous time scale 
is much longer than the infall time scale close to the black hole, the angular momentum 
remains almost constant, unless the viscosity is really very high (i.e. the a parameter 
close to unity). The supersonic region is well inside the accretion disk (typically within 50 
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Schwarzschild radii), far from the impact of any accretion stream that typically occurs at 
10 6 Schwarzschild radii. Therefore the supersonic inflow is expected to be approximately 
axisymmetric. Small non-axisymmetric perturbations may occur, of course, which is exactly 
the motivation for our study. 

The paper is structured as follows. Sections 2 and 3 describe the analytic and numerical 
models, respectively, the results of the numerical simulations are presented in Section 4, 
and conclusions with our plans for further investigations are discussed in Section 5. 



2. Modeling the problem 

We solve numerically the hydrodynamic equations describing the time dependent 
motion of a rotating, inviscid, ideal gas falling into a Schwarzschild black hole. We assume 
an initial axisymmetric configuration and a permanent axisymmetric boundary condition 
inflow. Since our study concerns the inner regions of the accretion disk (r < 50r g ) we may 
disregard the possible perturbation due to spiral shocks (Yukawa et al. 1997). Indeed, the 
validity of the extension of their result interior into the disk is far to be clear. Furthermore, 
the phenomenon, we discuss, may occur in accretion onto black holes in galactic nuclei, 
where a fully axisymmetric boundary condition may easily occur. To allow for analytical 
comparison and to be able to study the phenomenon with reasonable computer time, the 
motion is restricted to the z = plane and we assume that gradients in the z direction are 
negligible, i.e. d/dz = 0. The general relativistic effects are taken into account by using the 
Paczyhski & Wiita (1980) pseudo-Newtonian gravitational potential 

,, . GM c 2 
$(r) = - 



r 



9 



2(r/r g - 1) ' 



where r, G, c, M and r g are the radial distance, the gravitational constant, the speed of 
light, and the mass and the Schwarzschild radius (r g = 2GM/c 2 ) of the accreting object, 
respectively. This approximation is frequently used when accurate relativistic details are 
not required. The results show that both the initial axisymmetric and the final distorted 
shocks form at radial distances (r > 5r g ) for which the pseudo-Newtonian potential is fairly 
accurate. 

We solve the equations of hydrodynamics in polar coordinates r and <p for mass 
density p, radial momentum density pv r , angular momentum density rpv^ and energy 
(kinetic+thermal) density e. The four equations express conservation of mass, radial 
momentum balance, conservation of angular momentum, and energy balance, respectively. 

drp drpv r dpv v _ 

dt + dr + dip ~ U 1 } 
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dr 2 




2 d<S> 
pv^+p- rp— 







-rpv r — 
or 



d<5> 



(4) 



(5) 



(3) 



where the thermal pressure is p = (7 — l)[e — piyl + v^)/2] with the polytropic index 
7 = 4/3, which is appropriate for relativistic gas. We adimensionalize our equations by 
using the Schwarzschild radius and the speed of light as reference units, i.e. r g — 1 and 



The outer boundary of the computational domain is set far away from the BH at a 
radius r max (typically r max ~ 20 — 50r ff ), but closer than the outer sonic point, so that the 
inflow is supersonic. The inner boundary is chosen at a distance r m i n inside the inner sonic 
point where the flow behind the shock becomes supersonic again due to the gravitational 
acceleration. Typically we take r min ~ 1.5r g . The boundary conditions are periodic for the 
< ip < 2tt coordinate. 

Since the equations are linear in the mass and energy densities, we may choose the 
density of the infalling gas at the outer boundary of the computational domain to be unity, 
p(r max ) = 1. Once the density scale is fixed, the steady state solution is fully determined 
by three conserved quantities: the accretion rate A = rpv r , the specific angular momentum 
A = rv,p and the Bernoulli constant B = $ + (e + p)/p (sometimes referred to as 'energy'). 
To match the full solution, which contains an outer sonic point outside the domain, the 
accretion rate cannot be chosen independently of the other two constants, rather it should 
be calculated from the set of the algebraic equations expressing the mass conservation, the 
energy conservation and the polytropic equation of state following the procedure outlined 
in the Appendix. The remaining two parameters, the angular momentum and the Bernoulli 
constant can be chosen such that the steady shock is positioned at a distance r s h oc k- The 
dependence of the shock position on A and B is shown in Fig. |l| in the parameter range 
of our test cases. The supersonic inflow parameters are independent of the azimuthal 
coordinate if and, if not perturbed, the initial steady state flow would remain axi-symmetric 
forever. 

The initial condition is the steady solution perturbed a few units upstream of the shock 
at r perturb ~ r shock + 3r g and <p = 7r in a small region of size Ar = r perturb Aip = 0.05r perturb . 
We typically perturb the pressure or the density by a few per cent only, therefore the 
energy and momentum content of the perturbation is completely negligible relative to the 
respective quantities for the full domain. The perturbation is advected into the BH with 



c = 



1. 
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the flow as shown by the upper panels of Fig. |^. 

3. Numerical scheme 

To integrate numerically the hydrodynamic equations (0)-©, we use two different 
programs, VAC (Toth 1996, 1997) and another code implemented by Kuznetsov (Kuznetsov 
et. al. 1998). Both codes use Total Variation Diminishing (TVD, Harten 1983) type 
schemes. In VAC we used a characteristic based Lax-Wendroff type TVD scheme with a 
Roe type approximate Riemann solver (Roe 1986) and minmod or Woodward flux limiters 
(see Toth & Odstrcil 1996 for details), while in the other code a Lax-Friedrichs type 
scheme is applied with Osher-Chakravarthy anti-diffusion term (Chakravarthy & Osher 
1985). The combined Lax-Friedrichs-Osher scheme was suggested by Vyaznikov, Tishkin, 
& Favorsky(1989). Here the r and ip fluxes are added at the same time, while VAC uses a 
Strang type dimensional splitting. Both schemes are 2nd order accurate in space, the VAC 
code is temporally second order too, while the other code is only first order accurate in 
time. There are further differences in the grid settings and in the implementation of many 
numerical details. We have checked for several cases that the two codes give essentially 
identical results. 

The analytical equations are discretized in conservation form. In both codes, the 
physical and numerical fluxes are calculated for the angular momentum density rpv^ 
instead of the usual tangential momentum density pv v to achieve better angular momentum 
conservation. In VAC, this does not require the modification of the approximate Riemann 
solver itself, only the obtained fluxes are handled differently. 

Another concern is the temporally accurate evaluation of source terms on the right 
hand sides of the equations (§) and @, which arise from the curvature of the coordinate 
system and the gravitational force of the BH. The Lax-Friedrichs-Osher code is only first 
order accurate in time, thus very small time steps are used (Courant number « 0.1). In 
VAC, we follow the algorithm described by Ryu et. al. (1995) in their Appendix with a 
correction of their equation (A22), which should read q n+1 / 2 = (g n + q% dro + S n )/2, to 
achieve second order accuracy. This allows larger time steps with Courant numbers ~ 0.8. 
A naive implementation of a Lax-Wendroff type TVD scheme in polar coordinates leads to 
first order temporal accuracy and very inaccurate steady state shock positions for Courant 
numbers close to unity. 

Boundary conditions are easily implemented with two layers of ghost cells surrounding 
the physical domain. Supersonic inflow and outflow can be realized with fixed and 
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extrapolated ghost cell values, respectively. The fixed cell values are derived from the flow 
parameters: the accretion rate A, angular momentum A, and the Bernoulli constant B, 
which are listed for all the test cases in Table 1. Periodicity in the (p direction is trivial to 
maintain. 

Convergence of the numerical results was checked by redoing the same calculation with 
different grid resolutions N r x N v = 200 x 60, 200 x 120, 200 x 100, and 300 x 90. We 
also have the option of using a logarithmic grid spacing in the r direction (typically we 
use Arj + i/ Arj ps 1.003), which allows a better resolution of sharp gradients close to the 
BH. The results were found to be extremely similar on all grids. The evolution becomes 
qualitatively wrong only for a very coarse grid spacing 100 x 30, where the oscillations are 
suppressed by numerical diffusion and the perturbed shock returns to the symmetric steady 
state after a short transient. We may conclude that the results are not too sensitive to a 
small amount of dissipation. 

Although an analytical steady state can be easily derived by solving ordinary differential 
equations (Chakrabarti & Molteni 1993), it is better to obtain a numerical equilibrium, 
which is always slightly different from the analytical steady state solution due to the 
discretization errors. In particular, VAC can solve the time dependent axially symmetric 
one dimensional equations starting from a crude initial condition. The most efficient way 
to converge to steady state is to do a ID simulation with axial symmetry employing a fully 
implicit time stepping algorithm (see Toth, Keppens, & Botchev 1998 for details). 

The ID equilibrium solution is 'rotated' around the symmetry axis to get the 
starting configuration for the 2D runs, i.e. the 2D variables w^ 2D \i r ,i v ) = w^ lD \i r ) for 
w = p, pv-r^rpv^^e, where i r ,i v are the grid indices. The radial grid spacing must be the 
same for the 2D and ID grids, but the number of grid points in the tp direction can be 
chosen freely. We have tested that without the non-axisymmetric perturbation the 2D 
steady state obtained this way is stable and no fluctuations are introduced in the flow, since 
the numerical scheme used for the 2D simulation is the same as for the ID simulation. 

Once a steady solution with the shock is obtained, the small perturbation described 
at the end of the previous section is added in the cells around r = r pertur b, ip = n and the 
flow is evolved with the 2D code using explicit time stepping for hundred thousands of time 
steps. 
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4. Results 

Several test cases were studied in the parameter space denned by the specific angular 
momentum A = rv^ and the Bernoulli constant B = $ + (e + p)/p- Beside A and B, Table 
1 also lists the scaled accretion rate A/p(r max ) = —rv r p/p(r max ), which is necessary for 
setting up the boundary conditions, and the shock position r s h oc k for the axially symmetric 
steady state, which can be obtained by the procedure described in the Appendix. The cases 
are ordered by the radial shock distance and the specific angular momentum: the bigger 
case numbers (1, 2, ... 5) correspond to larger r s h oc k, and the subcases (a, b, ...) have larger 
and larger A, respectively. The outer boundary is located at r max = 56, 20, 25, 50 and 56 
for the five cases respectively. 

To check for time variations, we monitor the flow for the maximum shock radial 
distance, for the shock position at ip = 0, and for the mass and thermal energy contents of 
different sectors and annuli. We also analyzed snapshots at different times (see Fig. 0), an d 
visualized the simulation with animations. 

All cases, except the stable case 4, evolve rather similarly. As the perturbation 
crosses the hot subsonic postshock region, the shock becomes slightly distorted, and it 
oscillates with a small amplitude. After this transient, the phases of the small oscillations 
synchronize, and a very clear distorted shock develops with the end closed back as shown 
by Fig. j|. Further details of the final shock structure can be studied on Fig. §. 

The new axially asymmetric shock, with a dominant m = 1 mode, reaches a new "quasi 
steady state" with finite radial extent (in case 5b the simulation had to be stopped when 
the shock touched the outer boundary of the computational domain). The radial distance 
range of the twisted shock is not very different from the axially symmetric shock position 
r shock- The quasi steady state means that the flow changes periodically, without any sign 
of damping or instability. In case of "regular oscillations" the flow pattern rotates like a 
solid body with a period P, so it is a true steady state in the co-rotating coordinate system. 
In other cases more than one frequencies are present, which means that the solid body 
rotation is superimposed with a genuine oscillation of the flow pattern. When the ratio of 
the two periods are close to a small integer, "beating" can be observed in the monitored 
quantities. When many periods are present, the time variations become irregular, although 
the amplitude remains bounded. Table 1 summarizes the results and the observed periods 
for all cases. 

Certain trends can be easily identified from the table. For a fixed shock distance r s h oc /c, 
the larger the specific angular momentum A is, the more and usually longer periods are 
observed. This can also be seen in Fig. |5], where the time variation of the total mass in a 
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angular sector < ip < 2tt/3 is plotted. The figure contains further information. Note that 
the higher A is, the larger the average mass in the sector becomes, i.e. the further away 
the deformed shock is from the BH, since most of the mass is contained in the compressed 
post shock region. Furthermore, the amplitude of the mass variation also increases with A, 
which corresponds to a larger radial extent of the asymmetric shock pattern. Even in the 
most unstable case Id, however, the total mass in the sector increases by about 30% only. 

A clean example of the "beating" phenomenon can be seen in Fig. || which shows the 
total mass variation in the same angular sector for case 2c. The Fourier spectrum of the 
time variation, obtained with a discrete FFT, shows two well defined peaks with almost 
identical amplitudes in Fig. ^. Since the ratio of the frequencies is close to 2, strong beating 
can be observed in the time variation. 

The periods P are in the same range 60 < P < 280 for all cases listed in Table 1, 
despite the large variation in the rotation period at the shock distance, which is 2nr 2 shock /\ 
and it varies roughly from 100 to 2000 for cases 1 to 5, respectively. The radial extent of 
the deformed shock is not too large in most cases, e.g. for cases 2a and 3 the radial shock 
distance varies between 9 to 11 and 13 to 16, respectively. In case 5b, however, the shock 
has reached the outer boundary, at r max = 56, which is more than twice the shock distance 
t shock = 23.4 of the axially symmetric steady state. 

5. Conclusions 

We find that the axisymmetric shocks predicted by Chakrabarti's theory (Chakrabarti 
1990) in a rotating inviscid accretion flow are generally unstable to azimuthal perturbations, 
but the instability saturates at a low level, and a new, stable, asymmetric configuration 
develops with a strongly deformed shock rotating steadily. It is possible that shocks produced 
with very large angular momentum are unstable, in the sense that the perturbation triggers 
large deviations from the circular symmetry that extend up to the outer sonic point, but 
this does not happen in the test cases presented here. 

The new disk is no longer axisymmetric despite that the boundary conditions are the 
same as initially: axially symmetric supersonic inflow at the outer and supersonic outflow 
at the inner boundary. According to the numerical simulations, the asymmetric shock 
configuration is continuously self sustained and self reproducing around the BH. In general, 
this asymmetric configuration of the flow speed, density, and temperature produces time 
variations in any measured quantities. For a fixed shock distance we find that the time 
variations are regular periodic for low angular momentum flows and they are more irregular, 
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containing more frequencies, for flows with a larger angular momentum. 

The time variation could be observed and it could be a signature of the black hole's 
presence. As it can be seen in Table 1, typical periods are in the range 60 to 300 in the 
dimensionless units. Taking P = 100 as a typical case, we may convert it to physical units 
by multiplying with 2GM/c 3 ps 1O~ 5 M/M sec. For a black hole of ten solar mass, this 
gives approximately 0.01 seconds, while a BH with M = 1O 6 M would produce oscillations 
typically with 1000 second periods. 

We point out that although the general scenario seems different, there are similarities 
with the nonaxisymmetric disk instabilities studied by Blaes & Hawley (1988). Of course 
other important physical ingredients have to be included for a more realistic study of the 
shock behavior: full 3-dimensional treatment, physical viscosity, true cooling mechanism 
etc. Fully three dimensional simulations using the Versatile Advection Code of the same 
problem are in due course and will be presented in a forthcoming work. The role of the 
physical viscosity should also be further investigated, since in general, as shown in a similar 
context (cfr. Lanzafame, Molteni, & Chakrabarti 1998), but in axisymmetric conditions, 
moderate viscosity may produce changes in the shock structure and induce oscillations, 
while a large viscosity may stabilize the flow in a Keplerian regime, eliminating the shocked 
solutions. Our convergence studies indicate that the azimuthal instability is not very 
sensitive to a small amount of "numerical viscosity" . In any case we note that it has already 
been shown by time independent studies that the inner part of canonical accretion disks 
may have supersonic flows even for the viscosity parameter a about 0.01 (Chakrabarti, 
1996) or even larger (Igumenshev, Abramowicz, & Novikov 1998). 

Our investigation is intended as a first exploration of the new axially asymmetric 
solutions. Despite the many simplifying assumptions, we suggest that such a phenomenon 
may occur easily in real physical systems when viscosity is small and initial angular 
momentum is sub-Keplerian since the parameter range that leads to shocked solutions is 
fairly large and the perturbation required to trigger the asymmetric configuration is very 
small. The azimuthal instability described in this paper, in any case, is a fine example of 
how non linearity of the fluid dynamic equations may break the symmetry of the initial and 
boundary conditions. 
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A. Algebraic Equations for ID Steady State Solution 



Here we give a method to calculate the ID steady state solution, in particular the 
scaled accretion rate A/ p{r max ) and the shock position r s h oc k as a function of the specific 
angular momentum A and the Bernoulli constant B. This method requires (iterative) 
solutions of algebraic equations without the need to integrate ordinary differential equations 
(cfr. Chakrabarti 1989, 1990 for a method using integration). 

Let us introduce the sound speed a = ("fp/p) 1 ^ 2 and the Mach number in the radial 
direction M. = —v r /a. In case of an axially symmetric one dimensional steady state, i.e. 
d/dt = d/dip = d/dz = 0, we can integrate (0), (|), and to obtain three conserved 



quantities, A = rpv r , A 



rv v , and B 



( e + P)/p + $j respectively. For any continuous 



isentropic solution, the density and the sound speed are related as p = Ca 2 ^" / ~ 1 \ where C 
is a constant for the particular solution. We may now proceed to eliminate p, v r , v^, p, and 
e in favor of the conserved quantities A, B, C, and A, and the single unknown function, the 
Mach number Ai(r). First a 2 should be expressed from the Bernoulli equation 



B 



l,.., 2 A 2 



1 



- 2 (Ma) +_ + _- 



then it can be substituted into the mass conservation equation 

A 



-rpv r 



7+1 

CrMa-t- 1 



to arrive at the final equation for M(r) 



£ = f(M) ■ gx, B {r) 



with 



and 



f(M) 



g\A r ) = r 



M 



iM 2 



(7-1) 



7+1 
2(7-1) 



A 2 

B ~ 2T 2 + 277 



7 + 1 
2( 7 -l) 



(Al) 

(A2) 

(A3) 
(A4) 

(A5) 



The / function has a single maximum at M. = 1, and it can be inverted both in the 
subsonic < M. < 1 and supersonic M. > 1 regions. The g function has in general two 
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local minima at r% and r 2 with r\ > r 2 , which can be determined by solving the algebraic 
equation dg/dr = numerically. 

At large distances from the BH, the Mach number is M. <C 1, while close to the 
horizon M. 3> 1, thus any continuous solution has to have a sonic point with A4 — 1, where 
f[Ai(r)] has a maximum. Since the f(M) ■ g(r) product must be constant along the flow 
(|A3|), the maximum of / should be at one of the minima of g, i.e. at r\ or r 2 . Therefore we 
can have two isentropic solutions M.\ and Ai 2 



M 1>2 (r) = r 1 



(A6) 



9xA r ) 

where we use the subsonic branch of f~ l for r < r 12 and the supersonic branch for r > r 12 . 

The outer boundary conditions can now be easily determined from the Aii(r) solution, 
which connects the outer sonic point at r\ with the boundary located at r max . First 
•M-i^max) is calculated from the algebraic equation (|A6|) , then the sound speed a{r max ) 
from (|Al|), and finally the radial inflow speed v r = —Mia can be obtained. The scaled 
accretion rate is Aj p{r max ) = —r max v r (r max ). 



A standing shock can occur in the solution at r s h oc k if -A / li(r s / loc fc) 
M-2{r shock) < 1 are related by the Hugoniot relation 



> 1 and 



M 2 = h(Mi 



2 + ( 7 - l)M\ 



2 lV2 



2 1 m\ - (t - 1; 



(AT) 



First M.i{r shock ) can be determined by solving the algebraic equation 

/(■Mi) 9x, B (ri) 



(A8) 



next the shock position r shock = g x ^[/(l)5 f A,s( r 'i)/A / li(r s / locfc )] can be calculated. In general 
there can be two solutions, but only the outer one is stable (Nakayama 1994). 
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Fig. 1. — Shock location r shock versus the Bernoulli constant B = (e + p)/p + $ is plotted 
for different values of the specific angular momentum A = rv^. 

Fig. 2. — Contourlines of density are shown at four different times t. The simulation is done 
in r, ip coordinates, but the picture is shown in the x, y plane for easier interpretation. The 
initial perturbation has not reached the circular shock in the upper left snapshot, and it just 
crosses the shock front at t — 24. The new twisted shock appears after a long time t w 1000 
and it remains stable indefinitely. The parameters of this simulation are similar to those of 
case 2a in Table 1. 

Fig. 3. — Contour plot of the radial Mach number \v r \/a for the last snapshot (t = 13459) 
of Fig. [2| Contourlines in the subsonic region (Mach number below unity) are suppressed 
so that the shock location is shown as the last contour line of the incoming accretion flow. 
The inner region with unresolved contourlines corresponds to the supersonic region of the 
postshock infall. The Schwarzschild radius is indicated with the innermost dotted circle. 

Fig. 4. — Wire frame representation of pressure p and radial Mach number v r /a in the polar 
coordinates r, if. The black hole is at r = just outside the right boundary. For sake of 
clarity only 50 grid lines are drawn in both directions, but the actual resolution is 200 x 100. 
The shock happens to close back into itself at if « 7r. The density contourlines corresponding 
to this snapshot are shown in the left bottom panel of Fig. |2|. 

Fig. 5. — Time variation of the total mass in an angular sector < <p < 2tt/3 for cases la, 
lb, and Id, i.e. fixed shock distance r^o^ = 5.3 with increasing specific angular momentum 
A. The cases were run until a the amplitudes of the oscillations became bounded, thus the 
curve for case la is in the lower left part of the plot. 

Fig. 6. — Time variation of the total mass in an angular sector < <p < 2n/3. 

Fig. 7. — Fourier transform of the time variation shown in Fig. |6|. Note the two distinct 
peaks in the Fourier spectrum at P = 134 and P = 259 with almost equal amplitudes. The 
closeness of the period ratio to 1:2 is responsible for the "beating". 
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Table 1. Simulation Parameters and Results 



case 


T shock 


A/p(r max ) 


A 


B 


P 


comment 


la 


5.3 


4.0890 


1.7770 


.04333 


63 


regular osc. 


lb 


5.3 


4.4984 


1.8000 


.03220 


98 


regular osc. 


lc 


5.3 


4.7262 


1.8100 


.02701 


111,209 


beating 


Id 


5.3 


5.1041 


1.8225 


.02000 




irregular 


2a 


7.8 


3.0000 


1.8000 


.03630 


125 


regular osc. 


2b 


7.8 


3.0896 


1.8100 


.03170 


242,270 


beating 


2c 


7.8 


3.1834 


1.8200 


.02715 


134,259 


beating 


3 


12.7 


3.3450 


1.8200 


.03332 


215 


regular osc. 


4 


17.2 


4.2750 


1.8255 


.03331 




nearly stable 


5a 


23.4 


4.8766 


1.8620 


.02307 


273 


regular osc. 


5b 


23.4 


5.0442 


1.8720 


.02000 
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